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1. Introduction 

Over the past several years, intense experimental and theoretical interest has focused 
on the interaction of laser light with gaseous Bose-Einstein condensates (EEC's) PP. 
Examples include Bragg pulses applied to condensates and condensates confined in 
optical lattices|2j. Experimental and theoretical studies in these areas are important 
because condensates can be manipulated with Bragg pulses and strong optical lattices 
can confine a definite number of condensate atoms in each well. The Mott-insulator 
phase transition jSj can produce a system that might become a candidate for prototypical 
quantum computer. 

This paper describes a method for rapidly finding accurate solutions of the Gross- 
Pitaevskii equation describing a BEC in the presence of laser light. Such cases provide 
some of the most demanding numerical computations in the field of gaseous Bose- 
Einstein condensation. The "Hybrid" Lagrangian Method (HLM) described here is 
an instance of the Lagrangian Variational Method (LVM) a method for finding 
approximate solutions of the Gross-Pitaevskii (GP) equation, when the condensate 
wavefunction exhibits rapid variation along a single spatial direction. Examples of this 
include Bose-Einstein condensates in optical lattices and Bragg and standing-wave laser 
pulses to condensates. 

The method presented here is similar to several other methods presented in the 
literature jSj-jU]. The present method does have some significant differences. The most 
important difference is that the motion of the condensate transverse to the direction 
of the laser-light is coupled to the condensate motion along the light direction and 
vice versa. This coupling presents a challenge to find the correct, self-consistent, and 
variationally appropriate initial conditions. Such self-consistent initial conditions are 
not necessary in the other methods cited. 

This paper is organized as follows. The equations describing the basic Lagrangian 
Variational Method are derived by writing down the trial wavefunction and describing 
the meanings of the parameters. The steps required to arrive at the final HLM equations 
of motion are presented. The method is then applied to the case of a BEC in the presence 
of an accelerated optical lattice to show that it is possible to use a lattice to move a 
condensate to a different location in a non-adiabatic way. 

2. The Lagrangian Variational Method 

The standard LVM is a procedure for obtaining approximate solutions of the time- 
dependent Gross-Pitaevskii (GP) equation when a trial solution containing unknown 
parameters is assumed. This procedure can be summarized as follows. First, a 
Lagrangian density whose Euler-Lagrange equation is the GP equation is written down. 
Next, a trial wavefunction containing several time-dependent variational parameters is 
chosen. The Lagrangian whose generalized coordinates are these parameters is then 
obtained by inserting the trial wavefunction into the Lagrangian density and integrating. 
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The time-evolution equations for these parameters are derived using the usual Euler- 
Lagrange equations for the new generalized coordinates. 

It is possible to assume no particular form for the trial wavefunction and write down 
an Euler-Lagrange-type equation involving the Lagrangian density that yields the exact 
Gross-Pitaevskii partial differential equation. If one does assume a particular functional 
form for the trial wavefunction having time-dependent variational parameters, then 
equations of motion for these parameters are derived by integrating the Lagrangian 
density over all space yielding an ordinary Lagrangian depending only on the parameters 
and then applying standard Euler-Lagrange equations. 

The "Hybrid" Variational Method combines these two ideas in that the trial 
wavefunction assumed leaves the part of the wavefunction that depends on a particular 
coordinate completely unspecified while assuming a particular functional form (with 
time-dependent variational parameters) for the part of the wavefunction depending on 
the other coordinates. The resulting equation of motion for the unspecified part of the 
wavefunction is a partial differential equation while the equations of motion for the other 
variational parameters are ordinary differential equation that are first-order in time. 

The time-dependent GP equation for a confined condensate that is subjected to 
laser light has the following form. 

= -^^^^ + ^trap (r, t) ^ (r, t) + \4aser (r, t) ^ (r, t) + | ^ (r, t) I' ^ (r, t) . (1) 
The Lagrangian density, 
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L {r. d,r, dyr, o.r, Otr-, r, t) = l^h - ^*^) + 1;^ i^^i' 

+ (Krap (r,t) + V^aser (r,t)) |^|' + ^9 (2) 



(3) 



yields the GP equation when the variation 

= A ( \ — ( \ — ( \ — ( \ ^ 

Sip* dt\d{dtip*)J dx\d{dxip*)J dy \d (dyip*) J dz \d {dzip*) j dip* 

vanishes. This is the Euler-Lagrange equation that applies when the wavefunction is 
completely unspecified. 

To obtain a variationally approximate solution to the GP equation, we assume a 
trial wavefunction that contains a set of variational parameters. Each parameter is 
assumed to be a function of time. The most common example of a trial wavefunction 
is a three-dimensional gaussian: |i] 

iP {A (t) , w (t) , b (t) ; r) = A (t) e-^V«'^(t)+*/34i)^^e-yV«'?(t)+*/3«(t)?;2e-^V«'?(t)+i/3.(i)^^ (4) 

where w = {w^,Wy,w^) and b = {P^,Py,P^). 

By inserting this trial wavefunction into Eq. (jH)) and integrating over all space we 
obtain the ordinary Lagrangian: 

L{A{t),w{t),h{t))= J C[iP{A{t),w{t),h{t);r)]d\. (5) 
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In this way, we obtain a Lagrangian in which the variational parameters are the 
generahzed coordinates. Equations for the time dependence of these coordinates 
are obtained from the usual Euler-Lagrange equations. For example, for the width 
parameter Wx we would have the following. 

d_ fdL_\ _ ^ _ Q 

dt \dwxj dwx 

If one carries out this procedure for the gaussian trial wavefunction given above one 
finds, after a little algebra, the following equation of motion for Wx 
2 /i^ [W ah^N 

m'^Wx V / m^'WxWyWz 
Thus one obtains either a partial differential equation if the Euler-Lagrange equation 
for the Lagrangian density is used or an ordinary differential equation for a variational 
parameter if the integrated Lagrangian is used. The HVM combines these ideas and an 
example of the HVM will be derived next. 

3. The HVM equations of motion 

If a BEC is immersed in a one-dimensional optical lattice or subjected to optical 
laser pulses, then it is often the case that the condensate wavefunction exhibits spatial 
oscillations along the direction of propagation of the incident laser light that are much 
more rapid than along the transverse directions. In this case we can write a trial 
wavefunction with the following form: 

^ (0 (x, t) , wy (t) , wz {t) , (3y (t) , (t)) = <p [x, t) e-yyH(^)+^Mt)y\-^V2n.Ut)+iPAt>^ (g) 

Notice here that the variational parameters are 0, Wy, Wz, jSy, and jSz- We have chosen 
to represent the shape of the wavefunction in the transverse direction with a gaussian 
profile. In many cases the Thomas- Fermi approximation holds and the transverse density 
is more accurately represented by an inverted parabola. We shall see in section 13.3.21 
that the size of the variationally determined density agrees well with the Thomas-Fermi 
approximation when that limit applies. 

We shall assume that this trial wavefunction is normalized like this 

lijjf = N, (9) 

where N is the number of condensate atoms. The normalization condition constrains 
the values of the variational parameters 

/-f-OO P-\~C>0 /•-[-OO 

dx\(P{x,t)f dye'y'/< dze-'"'""' 
-oo J —oo J —oo 

+ 00 



ij^^'^Wy) (t^^'^Wz) / dx\(P{x,t)\^ = N. (10) 



Thus the wavefunction for the fast direction, (x, t), must satisfy 



,2 N 



dx\(t){x,t)\= . (11) 

-oo TTWyWz 

This can often be used to eliminate some variational parameters from the Lagrangian. 
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3.1. The Hybrid Lagrangian 

Since the functional dependence of the trial wavefunction on the y and z coordinates 
is assumed to be gaussian, we shall treat them in the standard LVM way. That is, we 
shall integrate the Lagrangian density only over these coordinates to get a "hybrid" 
Lagrangian 



L{(t)\d,<P\dt<P\Wy,w,,Py,p:)= dy j dzC[ilj{v,t)]. (12) 

J —oo 

where we take the full Lagrangian density to be 



1 / dip* n.dip\ h 



2 



|2 



c ir, d,r, dyr, o^r, dtr-, r, t) ^ -zn ^-^ip - j + ^ i^^i 

+ V,^,{v,t)\ipf + ^g\iP\* (13) 

and where 

Kxt (r, t) = Krap (r, t) + leaser (x, t) . (14) 

We will assume that the confining potential, Krap, has the form of an harmonic oscillator 
and that the potential created by the light field, Viaser, depends only on the x coordinate. 
The explicit forms of these are 

T/ /j-N ■'■ 22,-'- 22,-'- 22 

Vtrap(r,t) = -muj^x +-mujyy + -mu.z 

Maser (x, t) = - (1 + COS {2k^X - ^L^)) (15) 

where m is the atomic mass, [ujx^ujy^ujz) are the frequencies of the harmonic confining 
potential, /cl is the wavevector of the laser and 5l is the difference in frequencies of the 
two counterpropagating beams that comprise the laser pulse, VLq is the single-photon 
Rabi frequency of the condensate atom for the laser pulse, and A is the detuning of the 
incident laser pulse from atomic resonance. 

Calculating the form of the hybrid Lagrangian for the given trial wavefunction (Eq. 
(jHI) consists of (1) inserting the wavefunction into the full Lagrangian density, and (2) 
performing the integration over the y and z coordinates. We shall carry out these steps 
for each separate term of the full Lagrangian 

2m \ 2m J \ dx dx dy dy dz dz ) 

A = Kxt (r, t) \ip\^ = (^^mulx^ + ^mUyy^ + ^mulz'^ + Viaser {x, t)^ \ipf 

A = ^^?|^r (16) 

Then the resulting form of the hybrid Lagrangian for the given trial wavefunction is 
found from 

/+00 r+00 
dy / dzC,[ilj{r,t)]. 1 = 1,... A (17) 
-00 J —CO 
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where w = {wx^Wx) and j3 = {(3x,[3x)- The hybrid Lagrangian is found by inserting 
the trial wavefunction into the above and performing the required differentiations and 
integrations. 

The result is 




We now present the derivation of the equations of motion for the variational parameters. 



3.2. The Euler-Lagrange equations of motion 

The next step is to derive the equations of motion for 0, Wy, w^, Py, and Pz- We shall 
consider each in turn. We shall then discover that some further algebra will be required 
to make the resulting equations somewhat more convenient for calculation. 



3.2.1. The equation of motion for (j) (x, t) The Euler-Lagrange equation of motion for 
the variational parameter 0* is given by 



d I SLhybrid 



+ 



d dL 



hybrid 



dL 



hybrid 



dx 



0. 



(19) 



To evaluate this equation for our case we first compute the derivatives of ^hybrid with 
respect to 0*, 9f0*, and dxcf)* ■ Inserting these derivatives into Eq. (|Tn|) gives the following 
equation of motion for 0: 

h^ d\ 



ih— — I — ih 
dt 2 



Wy ^ Wz 
Wy Wz 



2m dx"^ 



laser 



(a;, t) + ^mulx"^^ (p + \(f)f (j) 



+ 



^ro 2 ^ / 1 

-npyW,, H — - 

2 ^ 2m I 2w?. 



+ ] + 7 



2„,,2 



niuo.w. 



y y 



1 



+ -h(3,wl + — , ^ 2 
2 2m \ 2wi 



By defining 



and 



1 



2 2 



x,y 



[/ext {X, t) = Vlaser {x, t) + ^mixllx^ 



(20) 



(21) 



(22) 
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we can write the equation of motion for as 



ih 



dt 



ih 



Wy 



+ 



+ i[fx + fy 



We can simphfy this equation considerably using the following transformation: 

(l){x,t)^N'/''4>{x,t) e-"W+*W 

where 

1 /•* 



The result is 



Wy (f) _^ {f 



Wy (f) W, (f) 

d4> 



dt' and b{t) = - [ {fy{t') + f,{t'))dt' 
Jo 



ih 



dt 



This equation can be simplified even further by noting that 

g-2a(t) ^ (0) W, (0) 
Wy (t) W^ (t) ■ 

The final form for the equation of motion for is the following 



A. , ^ 1 2 2\ 7 1 f^vi^y 



'w, (0)' 



2 ~ 



(23) 
(24) 
(25) 

(26) 
(27) 

(28) 



This is the equation of motion for 0. We now turn to the equations of motion for the 
widths and phases Wi and 



3.2.2. The equations of motion for the phases, (5i We turn now to the Euler-Lagrange 
equations for the gaussian phases. These equations will provide a relationship between 
the phases and the gaussian widths. These relationships will be useful in simplifying 
the equations for the widths. 

The Euler-Lagrangc (EL) equations for the gaussian phases, j3y and j3z, in terms 
of the hybrid Lagrangian arc as follows: 

Computing the required derivatives and inserting them into the EL equation gives 



(29) 



d_ 
di 
and 
d_ 
di 



■-hwl\<l^{x,t)f{ny'wy){ny'w,) 



hwl\cf>{x,t)f{'K'/'wy){ny'w,) 



h 



2 1 



2m 



2m 



(4/3,«;g|0(a;,i)|VV2^^,7rVV (30) 



(4/?,^2^|0(x,i)|VV2^,7rVV (31) 



We can simplify these equations by first considering the left-hand-side: 



d_ 
di 



■-hwl\cl>{x,t)\'(7r'/'wy){7T'/'wz) 



(32) 
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Thus we can write the EL equation for j3y as 



(4/5X) |0(x,t)|Vi/^.^,vri/V = |: (^K) {^'''^y) (^'^V) 



2m 



9t V2 

+ (^'^'^^) ('^'Ml (33) 

To simphfy this further, multiply both sides of the above equation by dx and integrate: 

-.2 



dx 



-1/2, 



TT ' Wy] ' W; 



-1/2, 



X 



+ 2K 



X 



dx ^ 
dt 



(34) 



We assume that the derivative and the integral can be interchanged in the last term on 
the right. Given this to be true, we can use the normalization condition (Eq. ()l()|l ) 

(35) 



(n'/^Wy) (n'/^w,) dx\(l){x,t)f = N 
to show that the last term in Eq. ()34|) is zero and thus we have 



(-) 



(-hwl 



dt V2 



(4A 



yWy 



2(3yWy 



d_ 
di 



(Wl) = 2WyWy 



(36) 



The EL equation for (3^ is similar and so we finally have 



m\ w. 



y 



2hJ w. 
m\ Wz 
2hJ Wz 



y 



(37) 



These equations can be used to eliminate Pi and Pi {i = x, y) in the expressions for the 



hh it) 



-npiW: H — 7; 

2 ^ ' 2m I 2w} 



1 



-mwiivi + 



^1 1 2 2 
2m2wf 4 * ' 



x,y. 



(38) 



This will be useful in simplifying the equations for the gaussian widths to which we now 
turn. 
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3.2.3. The equations of motion for the widths, Wi We next consider the EL equations 
for the widths Wi, where i = x,y. The EL equation for the Wi is given by 

d ( l9Lhybrid \ _ "^-^hybrid 

dt \ dwi ) dwi 



0. 



(39) 



We shall derive the equation for Wy] the derivation of the equation for Wz is similar. 
Inspection shows that ^hybrid is independent of Wy and so the above equation becomes 

f^i^hybrid 



dw,, 



0. 



(40) 



We can easily compute this derivative: 



dL 



hybrid 



ot \ \ 2m I ox ox 



+ 7 + l^^y< 



+ 



h$zwl\<i)f + 



2m ^ 
1 



+ V^aser ix,t) 
1 



— 

2m J \2wl 



|2 1 2 2 

I +^muj^x 



+ 2/3X^1 101^ + 



X 71 



-1/2 



1 

w: 



-mUyWy 



Defining 
Hx = ^ Im 



2m J dx dx 
we can rewrite Eq. as 

C^-^hybrid 



+ Viaser {x, t) 



+ 101^ + 7^101^ 



dw,, 



2 w,, \ 2m I 



1 



2wi ^ ^^^"^^ J ^ I'^^y'^y 



1 



2WyW^^ 



1 

-mu 

4 -w. 



2^1' 



y . 



H,,. 



Wy 

= 0. 

By multiplying by dx on both sides and integrating over x we obtain 



dL 



hybrid 



dw,, 



N{ -hPyWy + 



2m I 



2m J \2wyWz 



2p!^] + '-mul^' 

'WyJ 4 'Wy^ 



= 0. 



(41) 



(42) 



(43) 



(44) 



Hybrid Lagmngian Variational Method 
where 



10 



dx 



h Im 



dt 



+ 



' \ d(j)* d(t) 
2m J dx dx 



+ ^aser {x, t) 



(45) 



By using the equations of motion for the phases and their time derivatives we can 
ehminate the (3i and I3i from Eq. 



dw,. 



JM < -mwy + -m — Wz + 
4 ''Aw,, 



+ -mujyWy + -muj 



2^1 



2m J 

+ 7r{H,)wz = 



1 1 

+ 



2w'^ 2wywl^ 



(46) 



The equation for Wz is obtained by writing the above equation down with y and z 
interchanged. To simphfy these equations further we must consider the expression for 

We can express {Hx) in terms of 0, defined in Eq. fj24j) . by inserting that equation 
into Eq. (gHl) to get 

/4-CXD 2 
dx (47) 
-oo 



where 



dx Ihlm. 



dt 



+ 



' \ d(j)* 90 
2m I dx dx 



Vlaser {x, t) + -mUJ^X (f) + -qNc 



~ 2 



-2a{t) 



This can be significantly simphfied by using the equation of motion for 0. 



dx -gNe-'^^^'^ 
4 



Thus we can write a compact expression for (H^). Using Eq. (j47p we have 



hoo 1 , , ~ 4 

dx-gNe-^"^'^ 
4 



-hifyit) + fzit)) 



dx 



(48) 



(49) 



(50) 



Using the above expression we are at last ready to build the final equations of 
motion for the width parameters. The results are 



Wy + UyWy 

Wz + UjIwz 

where 



2\ 



,3 ~'~ ,,,2, 



m^ J w^Wz 



1 2A 
+ 



m^ / wi 



Wywf 



\U)^7^w'y{Q)wi{Q) 



• IgN 

dx 

4 m 



(51) 
(52) 

(53) 



This completes the full set of equations of motion for the gaussian parameters. 
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3.2.4- The final equations of motion Thus we may summarize the full set of equations 
of motion for this Hybrid Lagrangian Method. The trial wavefunction has the form 



ij (0, Wy, w,, Py, (5,) = <t> (x, t) e-^/V2«;g{^)+^/3.(^)y^e-.V2«;?(t)+^/?.(t)^^ 



(54) 



The variational parameters are (a;, t), the dependence of the wavefunction along the 
fast {x) direction, Wy and w^, the widths of the gaussians in the transverse dimensions, 
and (5y and (5^ the gaussian phases. The fast-direction wavefunction is written as 



^{x,t) 



^N^/^{x,t) 



X exp 



2 Jo 



Wy(f) Wj^ 
Wy{t') W,{t') 



dt' + lj\fy{t')+f,{t'))dt' 



(55) 



where 
hfi (t) 



2 ' 2m ylw]; 

1 



x,y. 



(56) 



and N is the number of condensate atoms. The equation of motion for is 



2m dx'^ 



+ Via.er {x, t) + -mwlx' (f) + -qN 



' Wy (0) \ ( Wz (0) ' 

Wy(t)) \W,{t)^ 



2 ~ 



The equations of motion for the widths are: 

2A 



Wy + UyWy 



1 

— + 



2A 

WyW^^'' 



where 



x{^)^.wiio)wuo)j_'yxi'-^ 



ct>. (57) 

(58) 
(59) 

(60) 



Once the widths are known, the phases can be obtained from the following relations: 



f m \ w. 



(61) 



Note that the equations for and the Wi form a closed set of equations that must be 
solved self-consistently. Now that we have derived the equations of motion, we next 
consider the method for numerical solution of these equations. 
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3.2.5. Scaled equations of motion To facilitate numerical solution of the HVM 
equations of motion it is important to cast them in scaled (dimensionless) units. The 
useful units of length {d) and time (l/o;) here are those of the geometrically averaged 
harmonic potential: 

1/2 



d = 



h 



1/3 



7i = 



UJi 



(62) 



Thus the x coordinate, gaussian widths, the time, and the wavefunction are scaled as 
follows 



X 



Pi f ■ \ 



ut 



4> 



(63) 



Transforming the HVM equations is straightforward. The results are as follows: the 
scaled equations of motion consist of the psuedo-lD GP equation, 



(rr / N 1 2 2 
f ( Viaser [x, t) + -^i^X 



dr 2 dx^ s 
the width equations. 



Attn 



' wy (o) \ fwAoy 



^'0, (64) 



-// 2- 1 2A(r) 

Wy + lyWy = + " ^ 

Wy WyWz 

1 2A (t) 
^ + 



WyWi 



where the scaled A parameter is 

A (r) = ^nwl (0) (0) 
and the scaled phase equations. 



r 

J—c 



dx 



(f) (x, r) 



w 

y 

2w,, 



and 



w' 



2wz 



(65) 



(66) 



(67) 



These are the scaled HVM equations. The solution of these equations cannot be found 
unless we know the initial conditions. We next address the problem of determining these 
initial conditions. 



3.3. Initial conditions for the HVM equations of motion 

The initial conditions for the HVM equations are assumed to represent a Bose- Einstein 
condensate statically held in a magnetic trap. The proper initial values for the fast- 
direction wavefunction (0(x, 0)), the transverse widths {wy{0) and u;^ (0)), and the 
transverse phases {j3y (0) and j3z (0)) must be stationary when propagated forward in 
time with the HVM equations. Determining the correct initial values for these quantities 
is the most difficult part of solving the HVM equations because the system of equations 
for (f) {x, 0), Wy (0) and Wz (0) are nonlinearly coupled to each other and must be solved 
self-consistently. These equations will be developed below and an approximate analytic 
solution for an important special case will be derived. 
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3.3.1. The equations defining the HVM initial conditions For a magnetically trapped 
condensate, the time evolution of the full wavefunction has the form 



(68) 



where p, is the scaled chemical potential. After making this transformation, the initial- 
condition equations are obtained by setting all of the remaining time derivatives in the 
HVM equations to zero: 



/i0 (x, 0) 



1^:1 2-2 

2 9x2 ^ 2^"^ 



0(x,O) 



0(x,O) (69) 



the parameter /i is to be determined. The equations for the widths at r = are 

1 2A (0) 



%Wy (0) 

(0) 



«;3 (0) ^ wl (0) (0) 
1 2A (0) 



(70) 



(0) Wy (0) (0) • 

Note that the parameter A (0) contains initial values of the widths and this motivates 
the following definition 



A-kN 



d 



dx 



0(x,O) =aowl{0)wUO) 



(71) 



The equations defining the initial widths then become 

1 



7i^y (0) 

(0) 



wl (0) 



2aQW^ (0) 
2aQWy (0) . 



(72) 



Equations (jUIH) and (f?^ together must be solved self-consistently to obtain the 
initial conditions for the HVM equations. Note that, while Eq. (j69p seems not to depend 
on the widths, this is not the case because the normalization of the full wavefunction 
must also be satisfied. In terms of the independent variables found in the HVM equations 
this becomes 

2 1 



dx 



0(x, 0) 



Txwy (O)w;^ (0)' 



(73) 



3.3.2. An approximate analytic solution For the case of a magnetic-trap potential that 
is axially symmetric along the fast direction and when the Thomas-Fermi limit is valid, 
it is possible to find an approximate analytical solution to the equations that define the 
proper variational initial conditions. While this solution is only approximate, it is useful 
because the case it covers is a common, realistic experimental situation. 

When the trap potential is cylindrically symmetric, then 7^ = 7^ = 7j^ and, by 
symmetry, it follows that 



Wy (0) = Wz (0) = W_L. 



(74) 



Hybrid Lagrangian Variational Method 



14 



In this case, the two Eqs. (f72|) become identicah 



1 



7_lW_l = — y + 2q;ow± (75) 



w 



As long as 2ao < this simple equation has the following solution 

Wy (0) = (0) =w^= (7I - 2ao) . (76) 

This is only the first step, however, since the factor ao is not known. 
The factor a^^ defined as 



"0 = 



, a 

A'kN ^ 



-00 _ 4 



dx 

00 



, 0(a;,O) , (77) 

evidently depends on the fast-direction wavefunction, (x, 0), which is, as yet, unknown. 
A simple, approximate expression for (x, 0) can be found, in the Thomas-Fermi limit, 
by neglecting the kinetic-energy term in Eq. ()69|) . 

By parameterizing the chemical potential in terms of the (also unknown) condensate 
radius along the fast direction 

= (78) 

and neglecting the kinetic energy in Eq. (jU^ . an expression for (x, 0) can be written 

as 

and zero otherwise. This formula can be used to obtain an expression for ao in terms 
of the condensate radius, xq\ 



(«)(f 



4-5 



Now both the initial widths and fast-direction wavefunction have been parameterized 
in terms of the condensate radius along the fast direction. 

The value of xq is determined by the normalization condition given in Eq. ()73|1 . 
Inserting the width and fast-direction wavefunction expressions in terms of Xq into the 
normalization condition gives an equation that determines this quantity: 

By solving this equation for xq, an approximate initial variational wavefunction can 
be obtained that will be stationary when propagated using the variational equations of 
motion derived earlier. 
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Figure 1. The average position of the condensate along the lattice direction. Here 
the lattice acceleration is 10 meters/ second^ and the turn-on times are Ti — 100 
and T2 — 800 /is (see text). 



4. Application: Moving a BEC with an optical lattice 

Finally, as an example of the usefulness of this method, it will be applied to the case 
of the behavior of a BEC in the presence of an accelerated optical lattice. It is possible 
to use an optical lattice to move a BEC from one place to another. The question arises 
as to whether a lattice can pick up a condensate, move it to a different location, and 
put it down at zero velocity. This is possible if the process is carried out completely 
adiabatically. However can it be accomplished rapidly? This question is considered in 
the following example. 

The GP equation will be solved for a condensate of = 10^ ^^Rb atoms confined 
in a cylindrically symmetric magnetic trap with frequencies u!±_ = {2tt) x 90 Hz and 
Uz = (27r) X 9 Hz. A stationary optical lattice is then turned on. The lattice lasers 
operate at a wavelength of A = 795 nm and its amplitude is ramped up from zero to a 
maximum of 5 ErecoU during a timespan Ti. A standing-wave lattice can be converted 
into a running-wave lattice by detuning one of the counterpropagating lasers from the 
other one. After the lattice is ramped up to its maximum value, the lattice is accelerated 
and allowed to run at a maximum velocity over a time T2 and then it is decelerated and 
ramped off symmetrically to the way it was turned on. The result of this sequence of 
events is illustrated in Fig. [H 

The position of the condensate is gauged by computing the average value of x 



_ J d^rj)* (r) xj) (r) _ /f^ dxcj)* (x) x0 {x) 
~ J d^r^j* (r) ip (r) ~ dx(j)* (x) (x) ' 

The second equality holds as a result of the assumed trial wavefunction. 
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The results of this calculation are shown in Fig. This figure shows that, while a 
stationary lattice is ramped on, the average position of the condensate is also stationary. 
When the lattice begins to move, the condensate is dragged along by fits and starts 
to a new position. These wiggles along the way are a result of the fact that the 
condensate is not a rigid body. When the lattice becomes stationary again its mean 
position moves with constant velocity but that this velocity is opposite the direction 
of the lattice motion. Further investigation shows that this velocity can be in either 
direction depending on Ti and T2 and that it is possible to use the lattice to pick up the 
condensate, move it to a new place, and put it down again with zero velocity, all done 
non-adiabatically. Although the condensate's final average velocity is zero, it will not 
be in its ground state and may be undergoing vigorous excitations. These results will 
be reported elsewhere. 

In conclusion this paper has presented a method for rapidly finding accurate 
approximate solutions of the GP equation for cases where a BEC is subjected to laser 
light. The 3+1 GP partial differential equation is reduced to a 1+1 partial differential 
equation plus a set of three second-order ordinary differential equations in time. The 
latter set of equations can be solved efficiently even on a grid that can accurately 
represent the solution along the rapidly varying laser-light direction. This method 
can be used in many applications involving condensate behavior under the influence of 
laser light. 
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